#include "CSpaceCharge.h"
#include "CBeam.h"
#include <iostream>
#include <fstream>
#include <cmath>
#include <stdio.h>
#include <iomanip>
#include <sstream>
#include <iterator>
#include <public.h>
#include "CVector3d.h"
using namespace std;
int xtimes=1;
int ytimes=1;
//#include "common.h"
//interpret() longitudinal ez  !!problem 1
#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])
//电流单位是mA,长度的单位是cm
static	CVector3d gCV;
CSpaceCharge::CSpaceCharge()
{
	//beam=new CBeam();
	cord=new vector<vector<double> >;                       
}
CSpaceCharge::~CSpaceCharge()
{
  //delete beam;
  delete cord;
  //delete cv;
}
void CSpaceCharge::formGrid()
{
	vvvdaa.resize(nx);
	vvvdex.resize(nx);
	vvvdey.resize(nx);
	vvvdez.resize(nx);
	compData.resize(nx);
	for(int i=0;i<vvvdaa.size();i++)
	{
		vvvdaa[i].resize(ny);
		vvvdex[i].resize(ny);
		vvvdey[i].resize(ny);
		vvvdez[i].resize(ny);
		compData[i].resize(ny);
		for(int j=0;j<vvvdaa[i].size();j++)
		{
			vvvdaa[i][j].resize(nz);
			vvvdex[i][j].resize(nz);
			vvvdey[i][j].resize(nz);
			vvvdez[i][j].resize(nz);
			compData[i][j].resize(2*nz);
			for(int k=0;k<vvvdaa[i][j].size();k++)
			{
				vvvdaa[i][j][k]=0;
				vvvdex[i][j][k]=0;
				vvvdey[i][j][k]=0;
				vvvdez[i][j][k]=0;
				compData[i][j][2*k]=0;
				compData[i][j][2*k+1]=0;

			}
		}  
	}

	vdxm.resize(nx-1);
	vdym.resize(ny-1);
	vdzm.resize(nz-1);
	vdxs.resize(nx);
	vdys.resize(ny);
	vdzs.resize(nz);
	xHalf=(nx-1)*dx/2.0;
	yHalf=(ny-1)*dy/2.0;
	zHalf=(nz-1)*dz/2.0;
	for(int i=0;i<nx-1;i++)
	{
		vdxs[i]=i*dx;
		vdxm[i]=vdxs[i]+dx/2.0;
	}
	vdxs[nx-1]=dx*(nx-1);
	for(int i=0;i<ny-1;i++)
	{
		vdys[i]=i*dy;
		vdym[i]=vdys[i]+dy/2.0;
	}
	vdys[ny-1]=dy*(ny-1);
	for(int i=0;i<nz-1;i++)
	{
		vdzs[i]=i*dz;
		vdzm[i]=vdzs[i]+dz/2.0;
	}
	vdzs[nz-1]=dz*(nz-1);
}
void CSpaceCharge::scheff()//only cal the space charge in grids. 
	//performing only when eht beam size changes a lot,for example the beam size becomes 5 times larger.
{
	double dXMax=beam->getXMax();
	double dYMax=beam->getYMax();
	if(dXMax>nx*dx/2-dZero||dYMax>ny*dy/2-dZero)
	{
		double dRMax=beam->getRmax();
		dx=dRMax*2*1.5/(nx-1);
		dy=dRMax*2*1.5/(nx-1);
		xHalf=(nx-1)*dx/2.0;
		yHalf=(ny-1)*dy/2.0;
		for(int i=0;i<nx-1;i++)
		{
			vdxs[i]=i*dx;
			vdxm[i]=vdxs[i]+dx/2.0;
		}
		vdxs[nx-1]=dx*(nx-1);
		for(int i=0;i<ny-1;i++)
		{
			vdys[i]=i*dy;
			vdym[i]=vdys[i]+dy/2.0;
		}
		vdys[ny-1]=dy*(ny-1);
	}
	cout<<"particle number is "<<iParticleNumTotal<<endl;
	if(dI==0) return;
	dQ =dI/(dFreq*iParticleNumTotal)*chargeOverMass;  // I is average current in one rf period

	double dxP,dyP,dzP,dphiP;
	int    ix1,ix2,iy1,iy2,iz1,iz2;
	double dx1,dx2,dy1,dy2,dz1,dz2;
	for(int np=0;np<cord->size();np++)
	{
		dxP=(*cord)[np][0]+xHalf; 
		dyP=(*cord)[np][1]+yHalf;
		dzP=(*cord)[np][2]+zHalf;
		ix1=int(dxP/dx);iy1=int(dyP/dy);iz1=int(dzP/dz);
		ix2=ix1+1;iy2=iy1+1;iz2=iz1+1;

		if(ix1<0) {ix1=0;ix2=0;}
		if(ix1>=nx-1){ix1=nx-1;ix2=nx-1;}
		if(iy1<0) {iy1=0;iy2=0;}
		if(iy1>=ny-1){iy1=ny-1;iy2=ny-1;}
		if(iz1<0) {iz1=0;iz2=0;}
		if(iz1>=nz-1){iz1=nz-1;iz2=nz-1;}

		dx2=1.0;dy2=1.0;dz2=1.0;
		if(ix1!=ix2)
		{
			dx2=(dxP-vdxs[ix1])/(vdxs[ix2]-vdxs[ix1]);
		}
		if(iy1!=iy2)
		{
			dy2=(dyP-vdys[iy1])/(vdys[iy2]-vdys[iy1]);
		}
		if(iz1!=iz2)
		{
			dz2=(dzP-vdzs[iz1])/(vdzs[iz2]-vdzs[iz1]);
		}
		dx1=1.0-dx2;
		dy1=1.0-dy2;
		dz1=1.0-dz2;
		vvvdaa[ix1][iy1][iz1]=vvvdaa[ix1][iy1][iz1]+dx1*dy1*dz1;
		vvvdaa[ix1][iy1][iz2]=vvvdaa[ix1][iy1][iz2]+dx1*dy1*dz2;
		vvvdaa[ix1][iy2][iz1]=vvvdaa[ix1][iy2][iz1]+dx1*dy2*dz1;
		vvvdaa[ix1][iy2][iz2]=vvvdaa[ix1][iy2][iz2]+dx1*dy2*dz2;
		vvvdaa[ix2][iy1][iz1]=vvvdaa[ix2][iy1][iz1]+dx2*dy1*dz1;
		vvvdaa[ix2][iy1][iz2]=vvvdaa[ix2][iy1][iz2]+dx2*dy1*dz2;
		vvvdaa[ix2][iy2][iz1]=vvvdaa[ix2][iy2][iz1]+dx2*dy2*dz1;
		vvvdaa[ix2][iy2][iz2]=vvvdaa[ix2][iy2][iz2]+dx2*dy2*dz2;
	}
	/*
		 ofstream of("vvvdaa.TXT");
		 cout<<"output charge in grid"<<endl;
		 double  charge=0;
		 for(int i=0;i<vvvdaa.size();i++)
		 for(int j=0;j<vvvdaa[i].size();j++)
		 for(int k=0;k<vvvdaa[i][j].size();k++)
		 {
		 of<<vdxs[i]<<" "<<vdys[j]<<" "<<vdzs[k]<<" "<<vvvdaa[i][j][k]<<endl;
		 charge+=vvvdaa[i][j][k];
		 }
		 cout<<charge<<endl;
	of.close();
		 */
	for(int i=0;i<vvvdaa.size();i++)
		for(int j=0;j<vvvdaa[i].size();j++)
			for(int k=0;k<vvvdaa[i][j].size();k++)
			{
				vvvdaa[i][j][k]=vvvdaa[i][j][k]*dQ/dx/dy/dz*1.0e6;
				compData[i][j][2*k]=vvvdaa[i][j][k];
				compData[i][j][2*k+1]=0;
			}
	//of<<vdxs[i]<<" "<<vvvdaa[i][ny/2][nz/2]<<endl;
	cout<<" begin fieldmap"<<endl;
	EFieldMap();
	//FieldError();
}

void CSpaceCharge::EFieldMap()      //求解每个场点网格点电场值
{
	cout<<"begin FFT3D"<<endl;
	ofstream of;     //test
	FFT3D(&compData);
	double dAngle;
	cout<<"begin potential"<<endl;
	Potential();
	cout<<"begin IFFT3D"<<endl;
	IFFT3D(&compData);
	/*
	//test
		 of.open("potential");
		 for(int i=0;i<compData.size();i++)
		 for(int j=0;j<compData[i].size();j++)
		 for(int k=0;k<compData[i][j].size();k=k+2)
		 {
		 of<<vdxs[i]-xHalf<<" "<<vdys[j]-yHalf<<" "<<vdzs[k/2]-zHalf<<" "<<compData[i][j][k]<<endl;
		 }
		 of.close();
	//end test
		 */
	cout<<"cal grid field"<<endl;
	double dBUL,dBUR,dBDL,dBDR,dFUL,dFUR,dFDL,dFDR;
	for(int i=0;i<vdxm.size();i++)
		for(int j=0;j<vdym.size();j++)
			for(int k=0;k<vdzm.size();k++)                           //  zyx
			{                                               
				dBUL=compData[i][j+1][2*k];              //  BUL--back-up-left---    (k--j+1--i)
				dBUR=compData[i+1][j+1][2*k];            //  BUR--back-up-right--    (k--j+1--i+1)
				dBDL=compData[i][j][2*k];                //  BDL  back-down-left     (k--j----i
				dBDR=compData[i+1][j][2*k];              //  BDR  back-down-right    (k--j----i+1
				dFUL=compData[i][j+1][2*(1+k)];          //  FUL  forward-up-left    (k+1-j+1-i
				dFUR=compData[i+1][j+1][2*(1+k)];        //  FUR  forward-up-right   (k+1-j+1-i+1
				dFDL=compData[i][j][2*(k+1)];            //  FDL  forward-down-left  (k+1-j---i
				dFDR=compData[i+1][j][2*(1+k)];          //  FDR  forward-down-right (k+1-j---i+1
				vvvdex[i][j][k]=-(dBUR+dBDR+dFUR+dFDR-dBUL-dBDL-dFUL-dFDL)/4.0/dx*100;
				vvvdey[i][j][k]=-(dBUL+dBUR+dFUL+dFUR-dBDL-dBDR-dFDL-dFDR)/4.0/dy*100;
				vvvdez[i][j][k]=-(dFUL+dFUR+dFDL+dFDR-dBUL-dBUR-dBDL-dBDR)/4.0/dz*100;
			}
			/*
		 cout<<"begin printing field in grids"<<endl;		
		 ofstream ofField("efield.txt");
		 for(int i=0;i<vvvdex.size()-1;i++)
		 for(int j=0;j<vvvdex[i].size()-1;j++)
		 for(int k=0;k<vvvdex[i][j].size()-1;k++)
		 ofField<<vdxm[i]-xHalf<<" "<<vdym[j]-yHalf<<" "<<vdzm[k]-zHalf<<" "<<vvvdex[i][j][k]
		 <<" "<<vvvdey[i][j][k]<<" "<<vvvdez[i][j][k]<<" "
		 <<sqrt((vdxm[i]-xHalf)*(vdxm[i]-xHalf)+(vdym[j]-yHalf)*(vdym[j]-yHalf)+(vdzm[k]-zHalf)*(vdzm[k]-zHalf))<<" "
		 <<sqrt(vvvdex[i][j][k]*vvvdex[i][j][k]+vvvdey[i][j][k]*vvvdey[i][j][k]+vvvdez[i][j][k]*vvvdez[i][j][k])<<" "<<endl;
		 ofField.close();
		 */
}
void CSpaceCharge::FieldError()
{
	//建立一个球坐标系，
	double dGuss=0;
	double dphi,dtheta,dR;//单位是deg
	int iphi,itheta;
	double dxRadius,dyRadius,dzRadius;
	int ix1,ix2,iy1,iy2,iz1,iz2;
	double dx1,dx2,dy1,dy2,dz1,dz2;
	double ex,ey,ez;
	iphi=180;
	itheta=90;
	dphi=2*pi/iphi;
	dtheta=pi/itheta;
	dR=0.2;
	ofstream error("error.txt",ios::app);
	for(int i=0;i<itheta;i++)
		for(int j=0;j<iphi;j++)
		{
			dxRadius=dR*sin(i*dtheta)*cos(j*dphi);
			dyRadius=dR*sin(i*dtheta)*sin(j*dphi);
			dzRadius=dR*cos(i*dtheta);
			ix1=int((dxRadius+xHalf)/dx);
			iy1=int((dyRadius+yHalf)/dy);
			iz1=int((dzRadius+zHalf)/dz);
			ix2=ix1+1;
			iy2=iy1+1;
			iz2=iz1+1;
			if(vdxm[ix1]>dxRadius&&ix1!=0&&ix1!=nx-1)
				ix2=ix1-1;
			else if((vdxm[ix1]>dxRadius&&ix1==0)||(vdxm[ix1]>dxRadius&&ix1==nx-1))
				ix2=ix1;
			if(vdym[iy1]>dyRadius&&iy1!=0&&iy1!=ny-1)
				iy2=iy1-1;
			else if((vdym[iy1]>dyRadius&&iy1==0)||(vdym[iy1]>dyRadius&&iy1==ny-1))
				iy2=iy1;
			if(vdzm[iz1]>dzRadius&&iz1!=0&&iz1!=nz-1)
				iz2=iz1-1;
			else if((vdzm[iz1]>dzRadius&&iz1==0)||(vdzm[iz1]>dzRadius&&iz1==nz-1))
				iz2=iz1;
			dx1=1.0;dy1=1.0;dz1=1.0;
			if(ix1!=ix2)
				dx1=(dxRadius+xHalf-vdxm[ix1])/(vdxm[ix1]-vdxm[ix2]);
			if(iy1!=iy2)
				dy1=(dyRadius+yHalf-vdym[iy1])/(vdym[iy1]-vdym[iy2]);
			if(iz1!=iz2)
				dz1=(dzRadius+zHalf-vdzm[iz1])/(vdzm[iz1]-vdzm[iz2]);
			dx2=1.0-dx1;dy2=1.0-dy1;dz2=1.0-dz1;
			ex=dz2*(dy2*(dx2*vvvdex[ix1][iy1][iz1]+dx1*vvvdex[ix2][iy1][iz1])
					+dy1*(dx2*vvvdex[ix1][iy2][iz1]+dx1*vvvdex[ix2][iy2][iz1]))
				+dz1*(dy2*(dx2*vvvdex[ix1][iy1][iz2]+dx1*vvvdex[ix2][iy1][iz2])
						+dy1*(dx2*vvvdex[ix1][iy2][iz2]+dx1*vvvdex[ix2][iy2][iz2]));
			ey=dz2*(dy2*(dx2*vvvdey[ix1][iy1][iz1]+dx1*vvvdey[ix2][iy1][iz1])
					+dy1*(dx2*vvvdey[ix1][iy2][iz1]+dx1*vvvdey[ix2][iy2][iz1]))
				+dz1*(dy2*(dx2*vvvdey[ix1][iy1][iz2]+dx1*vvvdey[ix2][iy1][iz2])
						+dy1*(dx2*vvvdey[ix1][iy2][iz2]+dx1*vvvdey[ix2][iy2][iz2]));
			ez=dz2*(dy2*(dx2*vvvdez[ix1][iy1][iz1]+dx1*vvvdez[ix2][iy1][iz1])
					+dy1*(dx2*vvvdez[ix1][iy2][iz1]+dx1*vvvdez[ix2][iy2][iz1]))
				+dz1*(dy2*(dx2*vvvdez[ix1][iy1][iz2]+dx1*vvvdez[ix2][iy1][iz2])
						+dy1*(dx2*vvvdez[ix1][iy2][iz2]+dx1*vvvdez[ix2][iy2][iz2]));
			dGuss+=dR*dR*1.0e-4*sin(i*dtheta)*dphi*dtheta*sqrt(ex*ex+ey*ey+ez*ez);
			//    error<<dGuss<<endl;
		}
	error<<dx<<" "<<dy<<" "<<dz<<" "<<nx<<" "<<ny<<" "<<nz<<" "<<(cord->size())
		<<" "<<dR<<" "<<dGuss<<" "<<1-dGuss/10.46<<endl;
	error.close();
}
void CSpaceCharge::FFT3D(vector<vector<vector<double> > > *vvvdComp)
{
	gsl_fft_complex_wavetable *compwave;
	gsl_fft_complex_workspace *compwork;
	compwork = gsl_fft_complex_workspace_alloc (nz);
	compwave = gsl_fft_complex_wavetable_alloc (nz);
	double aTmp[2*nz];
	//	ofstream of("comp-befor-FFT");
	//	for(int i=0;i<(*vvvdComp).size();i++)
	//	for(int j=0;j<(*vvvdComp)[i].size();j++)
	//	for(int k=0;k<(*vvvdComp)[i][j].size();k++)
	//	of<<vdxs[i]<<" "<<vdys[j]<<" "<<vdzs[k]<<" "<<(*vvvdComp)[i][j][k]<<endl;
	//	of.close();
	for(int i=0;i<nx;i++)
	{              
		for(int j=0;j<ny;j++)
		{
			for(int k=0;k<nz;k++)
			{
				REAL(aTmp,k)=(*vvvdComp)[i][j][2*k];
				IMAG(aTmp,k)=(*vvvdComp)[i][j][2*k+1];
			}
			gsl_fft_complex_forward(aTmp,1,nz,compwave,compwork);
			for(int k=0;k<nz;k++)
			{
				(*vvvdComp)[i][j][2*k]=REAL(aTmp,k);
				(*vvvdComp)[i][j][2*k+1]=IMAG(aTmp,k);
			}
		}
	}
	gsl_fft_complex_workspace_free(compwork);
	gsl_fft_complex_wavetable_free(compwave);
	//getVvv(&compData,"comp-FFT-nx",0,0,0,0,0);
	compwork=gsl_fft_complex_workspace_alloc(ytimes*ny);
	compwave=gsl_fft_complex_wavetable_alloc(ytimes*ny);
	double bTmp[2*ny*ytimes];
	for(int i=0;i<nz;i++)
	{
		for(int j=0;j<nx;j++)
		{           
			for(int k=0;k<ny;k++)
			{
				REAL(bTmp,k)=(*vvvdComp)[j][k][2*i];
				IMAG(bTmp,k)=(*vvvdComp)[j][k][2*i+1];
			}
			for(int k=2*ny;k<2*ny*ytimes;k++)
				bTmp[k]=0;
			gsl_fft_complex_forward(bTmp,1,ny*ytimes,compwave,compwork);
			for(int k=0;k<ny;k++)
			{
				(*vvvdComp)[j][k][2*i]=REAL(bTmp,ytimes*k);
				(*vvvdComp)[j][k][2*i+1]=IMAG(bTmp,ytimes*k);
			}
		}
	}
	gsl_fft_complex_workspace_free(compwork);
	gsl_fft_complex_wavetable_free(compwave);
	//getVvv(&compData,"comp-FFT-nx-ny",0,0,0,0,0);

	compwork=gsl_fft_complex_workspace_alloc(nx*xtimes);
	compwave=gsl_fft_complex_wavetable_alloc(nx*xtimes);
	double cTmp[2*nx*xtimes];
	for(int i=0;i<nz;i++)
	{
		for(int j=0;j<ny;j++)
		{           
			for(int k=0;k<nx;k++)
			{
				REAL(cTmp,k)=(*vvvdComp)[k][j][2*i];
				IMAG(cTmp,k)=(*vvvdComp)[k][j][2*i+1];
			}
			for(int k=2*nx;k<2*nx*xtimes;k++)
				cTmp[k]=0;
			gsl_fft_complex_forward(cTmp,1,xtimes*nx,compwave,compwork);
			for(int k=0;k<nx;k++)
			{
				(*vvvdComp)[k][j][2*i]=REAL(cTmp,xtimes*k);
				(*vvvdComp)[k][j][2*i+1]=IMAG(cTmp,xtimes*k);
			}
		}
	}
	gsl_fft_complex_workspace_free(compwork);
	gsl_fft_complex_wavetable_free(compwave);
	//getVvv(&compData,"comp-FFT-nx--ny-nz",0,0,0,0,0);
}
void CSpaceCharge::FFT1D(vector<vector<vector<double> > > *vvvdComp)
{
	gsl_fft_complex_wavetable *compwave;
	gsl_fft_complex_workspace *compwork;
	compwork = gsl_fft_complex_workspace_alloc (nz);
	compwave = gsl_fft_complex_wavetable_alloc (nz);
	double aTmp[2*nz];
	for(int i=0;i<nx;i++)
	{              
		for(int j=0;j<ny;j++)
		{
			for(int k=0;k<nz;k++)
			{
				REAL(aTmp,k)=(*vvvdComp)[i][j][2*k];
				IMAG(aTmp,k)=(*vvvdComp)[i][j][2*k+1];
			}
			gsl_fft_complex_forward(aTmp,1,nz,compwave,compwork);
			for(int k=0;k<nz;k++)
			{
				(*vvvdComp)[i][j][2*k]=REAL(aTmp,k);
				(*vvvdComp)[i][j][2*k+1]=IMAG(aTmp,k);
			}
		}
	}
	gsl_fft_complex_workspace_free(compwork);
	gsl_fft_complex_wavetable_free(compwave);
}
void CSpaceCharge::IFFT3D(vector<vector<vector<double> > >* vvvdComp)
{
	//getVvv(&compData,"comp-before-IFFT",0,0,0,0,0);
	gsl_fft_complex_wavetable *compwave;
	gsl_fft_complex_workspace *compwork;
	compwork=gsl_fft_complex_workspace_alloc(nx);
	compwave=gsl_fft_complex_wavetable_alloc(nx);
	//  double aTmp[2*nx*xtimes];
	double aTmp[2*nx];
	//  for(int i=0;i<2*nx*xtimes;i++)
	//  aTmp[i]=0;
	for(int j=0;j<ny;j++)
	{
		for(int i=0;i<nz;i++)
		{
			for(int k=0;k<nx;k++)
			{
				//        REAL(aTmp,xtimes*k)=(*vvvdComp)[k][j][2*i];
				//        IMAG(aTmp,xtimes*k)=(*vvvdComp)[k][j][2*i+1];
				REAL(aTmp,k)=(*vvvdComp)[k][j][2*i];
				IMAG(aTmp,k)=(*vvvdComp)[k][j][2*i+1];
			}
			gsl_fft_complex_inverse(aTmp,1,nx,compwave,compwork);
			for(int k=0;k<nx;k++)
			{
				(*vvvdComp)[k][j][2*i]=xtimes*REAL(aTmp,k);
				(*vvvdComp)[k][j][2*i+1]=xtimes*IMAG(aTmp,k);
			}
		}
	}
	gsl_fft_complex_workspace_free(compwork);
	gsl_fft_complex_wavetable_free(compwave);
	//getVvv(&compData,"comp-IFFT-nx",0,0,0,0,0);
	compwork = gsl_fft_complex_workspace_alloc (ny);
	compwave = gsl_fft_complex_wavetable_alloc (ny);
	double bTmp[2*ny];
	//memset(bTmp,0,2*ny*ytimes*sizeof(double));
	for(int i=0;i<2*ny;i++)
		bTmp[i]=0;
	for(int i=0;i<nx;i++)
	{
		for(int j=0;j<nz;j++)
		{
			for(int k=0;k<ny;k=k++)
			{
				REAL(bTmp,k)=(*vvvdComp)[i][k][2*j];
				IMAG(bTmp,k)=(*vvvdComp)[i][k][2*j+1];
			}
			gsl_fft_complex_inverse(bTmp,1,ny,compwave,compwork);
			for(int k=0;k<ny;k++)
			{
				(*vvvdComp)[i][k][2*j]=ytimes*REAL(bTmp,k);
				(*vvvdComp)[i][k][2*j+1]=ytimes*IMAG(bTmp,k);
			}
		}
	}
	gsl_fft_complex_workspace_free(compwork);
	gsl_fft_complex_wavetable_free(compwave);
	//getVvv(&compData,"comp-IFFT-nx-ny",0,0,0,0,0);

	compwork = gsl_fft_complex_workspace_alloc (nz);
	compwave = gsl_fft_complex_wavetable_alloc (nz);
	double cTmp[2*nz];
	for(int i=0;i<nx;i++)
	{
		for(int j=0;j<ny;j++)
		{
			for(int k=0;k<nz;k=k++)
			{
				REAL(cTmp,k)=(*vvvdComp)[i][j][2*k];
				IMAG(cTmp,k)=(*vvvdComp)[i][j][2*k+1];
			}
			gsl_fft_complex_inverse(cTmp,1,nz,compwave,compwork);
			for(int k=0;k<nz;k++)
			{
				(*vvvdComp)[i][j][2*k]=REAL(cTmp,k);
				(*vvvdComp)[i][j][2*k+1]=IMAG(cTmp,k);
			}
		}
	}
	gsl_fft_complex_workspace_free(compwork);
	gsl_fft_complex_wavetable_free(compwave);
	//getVvv(&compData,"comp-IFFT-nx-ny-nz",0,0,0,0,0);
}
void CSpaceCharge::IFFT1D(vector<vector<vector<double> > >* vvvdComp)
{
	gsl_fft_complex_wavetable *compwave;
	gsl_fft_complex_workspace *compwork;
	compwork = gsl_fft_complex_workspace_alloc (nz);
	compwave = gsl_fft_complex_wavetable_alloc (nz);
	double cTmp[2*nz];
	for(int j=0;j<ny;j++)
	{
		for(int i=0;i<nx;i++)
		{
			for(int k=0;k<nz;k=k++)
			{
				REAL(cTmp,k)=(*vvvdComp)[i][j][2*k];
				IMAG(cTmp,k)=(*vvvdComp)[i][j][2*k+1];
			}
			gsl_fft_complex_inverse(cTmp,1,nz,compwave,compwork);
			for(int k=0;k<nz;k++)
			{
				(*vvvdComp)[i][j][2*k]=REAL(cTmp,k);
				(*vvvdComp)[i][j][2*k+1]=IMAG(cTmp,k);
			}
		}
	}
	gsl_fft_complex_workspace_free(compwork);
}
void CSpaceCharge::Potential()       
{     
	double dkx[nx],dky[ny],dkz[nz],dif3K[nx][ny][nz];
	double dTmpx,dTmpy,dTmpz;
	double dif;
	dkx[0]=0;dky[0]=0;dkz[0]=0;

	//FFT3D======================================
	double dLz=dz*(nz);
	double dLx=dx*(nx);
	double dLy=dy*(ny);

	//begin correct
	for(int k=0;k<ny;k++)
	{
		if(k<ny/2)
			dky[k]=k*(2*pi/dLy);
		else dky[k]=(k-ny)*(2*pi/dLy);
	}
	for(int i=0;i<nx;i++)           //when i=0.............!!!!!!!!
	{
		if(i<nx/2)
			dkx[i]=i*(2*pi/dLx);
		else dkx[i]=(i-nx)*(2*pi/dLx);
	}
	for(int j=0;j<nz;j++)
	{
		if(j<nz/2)
			dkz[j]=j*(2*pi/dLz);
		else       dkz[j]=(j-nz)*(2*pi/dLz);
	}
	//=====orgion====================
	for(int k=0;k<ny;k++)
		for(int i=0;i<nx;i++)           //when i=0.............!!!!!!!!
			for(int j=0;j<nz;j++)
			{

				if((i==0)&&(j==0)&&(k==0))
				{
					compData[i][k][2*j]=0;
					compData[i][k][2*j+1]=0;
					}
				else{
					if(k==0) dTmpy=1;
					else dTmpy=sin(dky[k]*dy/2.0)/(dky[k]*dy/2.0);
					if(i==0) dTmpx=1;
					else dTmpx=sin(dkx[i]*dx/2.0)/(dkx[i]*dx/2.0);
					if(j==0) dTmpz=1;
					else dTmpz=sin(dkz[j]*dz/2.0)/(dkz[j]*dz/2.0);
					dif=dTmpx*dTmpx*dkx[i]*dkx[i]*1.0e4+dTmpz*dTmpz*dkz[j]*dkz[j]*1.0e4+dTmpy*dTmpy*dky[k]*dky[k]*1.0e4;
					compData[i][k][2*j]=compData[i][k][2*j]/dif/dPermit;
					compData[i][k][2*j+1]=compData[i][k][2*j+1]/dif/dPermit;
				}
			}
}
void CSpaceCharge::setFreq(double _dFreq)
{
	dFreq=_dFreq*1e6;
	dWaveLength=clight/dFreq;
}
void CSpaceCharge::InputBeam()
{
	iParticleNumTotal=beam->getGoodParticle();
	cord->resize(iParticleNumTotal);
	double    dCenterZ=beam->getZCenter();
	double    dGammaAve=beam->getEnergyAve()/dErest+1.0;
	double    dXAve=beam->getAveSize()->getX();
	double    dYAve=beam->getAveSize()->getY();
	int       lossflag=0;
	double    dZ0;
	vector<double> vd(4,0);
	for(int i=1;i<cord->size()+1;i++)
	{
		lossflag=beam->p.at(i).getLossFlag();
		if(lossflag==0)
		{         
			dZ0=beam->p.at(0).getVelocity().getZ()*(beam->p.at(i).getT0()-beam->p.at(0).getT0());
			//(*cord)[i].resize(4);
			vd[0]=(beam->p[i].getPosition().getX())*100;                    //x    in cm
			vd[1]=(beam->p[i].getPosition().getY())*100;                    //y    in cm
			vd[2]=(beam->p[i].getPosition().getZ()-dCenterZ-dZ0)*100*dGammaAve;//dzcenter can not bo omited
			(*cord)[i-1]=vd;
			//it maters nothing wehther +dZ0 or -dZ0
		}
	}
	dI=beam->getCurrent()/1000.0;      
}
void CSpaceCharge::InputBeam(vector<vector<double> >* _vvdBeam,double _dI)
{
	cord=_vvdBeam;                                                    
	dI=_dI/1000.0;
	iParticleNumTotal=_vvdBeam->size()-1;      //第一个粒子是同步粒子，不包括在粒子数范围内
}
void CSpaceCharge::InGrid(int _nx,int _ny,int _nz,double _dx,double _dy,double _dz)
{
	nx=_nx*2;           //in cm
	dx=_dx;
	ny=_ny*2;
	dy=_dy;
	nz=_nz*2;
	dz=_dz;
}
CVector3d* CSpaceCharge::Interpret(double _dxP,double _dyP,double _dzP)  //输入单位是m
{
	double dxP,dyP,dzP;
	dxP=_dxP*100+xHalf;dyP=_dyP*100+yHalf;dzP=_dzP*100+zHalf;
	double dx1,dy1,dz1,dx2,dy2,dz2;
	int ix1,iy1,iz1,ix2,iy2,iz2;
	double ex,ey,ez;
	//	dxP=(*cord)[i][0]+xHalf;
	//	dyP=(*cord)[i][1]+yHalf;
	//	dzP=(*cord)[i][2]+zHalf;
	ix1=int(dxP/dx);
	iy1=int(dyP/dy);
	iz1=int(dzP/dz);
	ix2=ix1+1;
	iy2=iy1+1;
	iz2=iz1+1;
	if(vdxm[ix1]>dxP&&ix1!=0)
		ix2=ix1-1;
	else if(vdxm[ix1]>dxP&&ix1==0)
		ix2=ix1;
	if(vdym[iy1]>dyP&&iy1!=0)
		iy2=iy1-1;
	else if(vdym[iy1]>dyP&&iy1==0)
		iy2=iy1;
	if(vdzm[iz1]>dzP&&iz1!=0)
		iz2=iz1-1;
	else if(vdzm[iz1]>dzP&&iz1==0)
		iz2=iz1;
	dx2=1.0;dy2=1.0;dz2=1.0;
	if(ix1!=ix2)
		dx2=(dxP-vdxm[ix1])/(vdxm[ix2]-vdxm[ix1]);
	if(iy1!=iy2)
		dy2=(dyP-vdym[iy1])/(vdym[iy2]-vdym[iy1]);
	if(iz1!=iz2)
		dz2=(dzP-vdzm[iz1])/(vdzm[iz2]-vdzm[iz1]);
	dx1=1.0-dx2;dy1=1.0-dy2;dz1=1.0-dz2;
	ex=dz1*(dy1*(dx1*vvvdex[ix1][iy1][iz1]+dx2*vvvdex[ix2][iy1][iz1])
			+dy2*(dx1*vvvdex[ix1][iy2][iz1]+dx2*vvvdex[ix2][iy2][iz1]))
		+dz2*(dy1*(dx1*vvvdex[ix1][iy1][iz2]+dx2*vvvdex[ix2][iy1][iz2])
				+dy2*(dx1*vvvdex[ix1][iy2][iz2]+dx2*vvvdex[ix2][iy2][iz2]));
	ey=dz1*(dy1*(dx1*vvvdey[ix1][iy1][iz1]+dx2*vvvdey[ix2][iy1][iz1])
			+dy2*(dx1*vvvdey[ix1][iy2][iz1]+dx2*vvvdey[ix2][iy2][iz1]))
		+dz2*(dy1*(dx1*vvvdey[ix1][iy1][iz2]+dx2*vvvdey[ix2][iy1][iz2])
				+dy2*(dx1*vvvdey[ix1][iy2][iz2]+dx2*vvvdey[ix2][iy2][iz2]));
	ez=dz1*(dy1*(dx1*vvvdez[ix1][iy1][iz1]+dx2*vvvdez[ix2][iy1][iz1])
			+dy2*(dx1*vvvdez[ix1][iy2][iz1]+dx2*vvvdez[ix2][iy2][iz1]))
		+dz2*(dy1*(dx1*vvvdez[ix1][iy1][iz2]+dx2*vvvdez[ix2][iy1][iz2])
				+dy2*(dx1*vvvdez[ix1][iy2][iz2]+dx2*vvvdez[ix2][iy2][iz2]));
	//cv=new CVector3d(ex,ey,ez);
	//return cv; 
	gCV.setX(ex);
	gCV.setY(ey);
	gCV.setZ(ez);
	return &gCV;
}
